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<D Abstract 

43 "We introduce a fully generalized quiescent chemical reactor system in arbitrary space 

dim = 1,2 or 3, with n 6 N chemical constituents ctj, where the character of the numerical 
solution is strongly determined by the relative scaling between the local reactivity of species 
cti and the local functional diffusivity @ij{a) of the reaction mixture. We develop an operator 
time-splitting predictor multi-corrector RK-LDG scheme, and utilize /ip-adaptivity relying 
only on the entropy J?^ of the reactive system 9\. This condition preserves these bounded 
Qh nonlinear entropy functionals as a necessarily enforced stability condition on the coupled 

system. We apply this scheme to a number of application problems in chemical kinetics; 
I including a difficult classical problem arising in nonequilibrium thermodynamics known as 

the Belousov-Zhabotinskii reaction where we utilize a concentration-dependent diffusivity 
tensor 2$ij(a), in addition to solving a simple equilibrium problem in order to evaluate the 
numerical error behavior. 

in 
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§1 Introduction 

Chemical reactors are of fundamental importance in a large array of scientific fields, spanning 
applications in chemistry and chemical engineering [321 ESI ES] , mechanical and aerospace engi- 
neering [6T] , atmospheric and oceanic sciences [64] , astronomy and plasma physics [2S1 ES] ; as well 
as generally in any numbers of biologically related fields (viz. [55] for example). More fundamen- 
tally it is the basic prevalence of these dynamic reactive chemical systems in nature that makes 
the ability to effectively model them so essential. 

From a theoretical point of view, much of the underlying theory for reactor systems may be 
found in [15] 125] 135] , where generally reactor systems may be derived using kinetic theory by 
way of a Chapman-Enskog or Hilbert type perturbative expansion, which immediately raises a 
set of important concerns that are far beyond the present scope of this paper (see for example 
[68] for an example of the formal complications that may arise in rigorous treatments). Here we 
rather restrict ourselves to the study of a set of simplifications leading to a generalized system of 
reaction-diffusion equations, that may be referred to collectively as quiescent reactors. 

The foundational theory provides that quiescent reactor systems may be derived explicitly from 
fluid particle systems (i.e. the Boltzmann equation), where the reaction diffusion multivariate 
master equation [121 EHl EDJ serves as the rigorous justification underlying the model. From the 
point of view of the experimental sciences, a quiescent reactor may be defined as a chemical reactor 
system where no explicit stirring is either present or plays a significant role in the dynamic behavior 
of the medium. It must be noted however, that often in these quiescent experimental systems, 
as seen in [53], heat gradients may be utilized for chemical catalysis. Because of the complicated 
formal coupling between state variables, in particular those of density p = p(t,x), temperature 
i? = ■#(£, x) and internal energy $ = $(t,x), all of which may have the effect of imparting local 
velocity gradients on fluids elements, we restrict ourselves in this paper to isothermal systems, and 
define a quiescent reactor as one which, up to the transport properties of the system, is diffusion 
dominated in the sense of the Fickian regime (which we expand upon below); or, more precisely, 
may be modeled up to an implicit stochasticity as a system of reaction-diffusion equations [31 ES] ■ 

More clearly, from the point of view of the simplified mathematics of the system, such a 



restriction to the quiescent regime may be presented as any flow system obeying (2.1) which 
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satisfies the following approximate bound: 

V x -(a i u)<^' i \V x a i \ 2 + ^ i A x a i , (1.1) 
where u = u(t, x) is the flow velocity of the system, and where when the formal inequality holds 



(i.e. for < in 1.1), the quiescent approximation is particularly strong. Formally we can say that if 
the concentration and velocity gradients are comparable V x Uj ~ V x aj for each constituent i and 
each component j, or even more strongly whenever V x Uj < V x «j, then if the velocity components 



are bounded from above by Uj < a^'Vzlnaj — then (1.1) is satisfied, and strictly satisfied 
when the bounds are precise (i.e. < ==>- <). In the case of, for example, the Chapman-Enskog 
expansion of ^ as developed in §4, this merely suggests that for a bounded concentration gradient 
IV^c^I < C the diffusive gradient is controlled from below such that ^ > ko^ for k = n(aj^i, C) 
having only functional dependencies on the fractional weighting of the other constituents of the 
fluid. By contrast, when the diffusion coefficient is taken to be constant such that @- = 0, it 
follows that in local areas of appreciable concentration, i. e. oti 3> 0, the advection must scale with 
diffusive collisions, and similarly in areas of measurable velocities it is the concentration gradient 
which must scale with the collisional motions. 

From the point of view of the physics and chemistry of the system, a diffusion dominated flow 
regime is merely one in which the random collisional molecular motion of the fluid dominates 
the advective flow characteristic. Such systems are frequently used as approximate models to 
restrict to systems that implicitly contain substantially more complicated dynamics (e.g. such as 
in combustion models [72]). 

In fact, it is remarkable the number of complicated and important physical phenomena that are 
understood merely by way of modeling coupled reaction-diffusion equations. For example, the spa- 
tially distributed FitzHugh-Nagumo model is a reaction-diffusion system of primary importance 
in tracking the formation, propagation and recovery of action potentials in biological and artificial 
neural networks [2T] • In fact the Nagumo formulation [51] (of which the FitzHugh-Nagumo model 
may viewed as a special case) for single component reaction-diffusion models comprises the core 
of the underlying mathematics responsible for the chemical basis of morphogenesis in biological 
processes (e.g. the Kolmogorov-Petrovsky-Piskunov (KPP) Equation) [5"4"| 169]. Fisher's equation 
is also a biologically relevant reaction-diffusion system, and is used for modeling the propagation 
of genetic variation over sample populations [56] . while in plasma physics the modeling of multi- 
component reactive hot plasmas are often reduced down to systems of coupled reaction-diffusion 
equations [71] . Additionally, reaction-diffusion models are of central importance in the study of 
phase-field models j2j [66] and nonequilibrium thermodynamics (we provide a detailed discussion 
of the latter in §4.3), just to mention a very sparse few. 

One important and emergent feature of coupled reaction-diffusion systems is the nuance that 
arises in understanding that a diffusion dominated regime is not necessarily a diffusion limited 
regime - - e.g. in the sense of the standard parlance of analytic chemistry [70]. That is, the 
diffusion rates (or diffusion "velocity") will limit the reaction front in a reaction whose kinetics 
occur on shorter timescales than the particles diffuse (which is a diffusion limited process), but in 
a reaction with timescales that are appreciably longer than the timescales of the diffusion rates 
of the systems components, the chemical reaction rate can becomes the limiting process (i.e. a 
reaction limited process), and the diffusion regime switches from a diffusion limited process to a 
fully diffusion dominated process. In contrast, when the reaction rates occur much faster than 
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the diffusion rates, and when the domain is for example homogenized, then we see a fully reaction 
dominated process. Notice that different parts of the domain may be characterizes by different 
regimes. 

It should further be noted, as discussed for example in [19] and jl3], that frequently one 
encounters a full decoupling between the reactive and elastic regimes when the reactive time scales 
are much slower than those of the dissipative time scales. This does not, however, account for the 
popular engineering trend towards multiscale applications [9j [37] , where substantial differences in 
reaction and diffusing timescales may be present and yet still coupled through a standard reactor 
regime. 

In fact, it is precisely this difference in relative timescales between the diffusion processes of 
the system and the reaction processes of the system which makes developing a general numerical 
scheme difficult. One generally finds, when solving a parabolic system, that the stability condition 
on the timestep At < CAx 2 makes formulating an explicit solution unattractive, and implicit 
methods are favoured. However, in generalized quiescent reactor systems the timescales of the 
reactive components of the system may vary wildly in the effective scaling (viz. reactions rates 
on the order of 10 -15 to 10 15 in standard units), while the diffusion scaling may demonstrate 
substantially less variation. Because of this, the time-stepping limitations in quiescent reactor 
models may be either strongly reaction limited, or strongly diffusion limited, or both — in the 
sense of oscillating between the regimes, or being split across the regimes. That is, in reactor 
systems, where many different reactions may be occurring simultaneously, the depletion of a 
certain constituent a, at time t n may cause a local in time transition from a reaction dominated 
time-stepping regime to a diffusion dominated or diffusion limited time-stepping regime, and vice 
verse (in the sense of the operator-splitting regimes of [62J). Or, as may occur, parts of the domain 
may be dominated by constituents which are inert with respect to each other, are uncoupled, and 
which consequently operate with respect to fundamentally different regimes (i.e. diffusion versus 
reaction limited regimes). 

It must be additionally noted here that this distinction is in many ways a simplification of what 
can be a very subtle interplay between the reactive and diffusive modes present in complicated 
reactive mixtures. For example, it is well known that Fick's law of diffusion is in some cases ill- 
suited for describing the behavior of some hysteretic mixtures, or diffusion regimes with memory 
(e.g. such as the electrochemistry induced near an active electrolytic cell [3j [23], [26] ) . In these cases, 
the local propagation speed caused by the gradient of the concentration forces the Fick's component 
of the diffusion to obey a telegraph equation, which may often reduce to an integro-differential 
equation over all time [0, T] coupled to a mass transport equation in the reactive components. 
These complications arise in systems that demonstrate large variations in concentrations over 
short time frames, though a large class of reactions demonstrate even more complicated and 
subtle behavior that might require the inclusion of quantum effects, such as in |1S] . As a general 
rule we will not directly address these complications below, as we uniformly make the assumption 
that the reaction-diffusion system of equations employed is an appropriate approximate model for 
the system in question. 

Nevertheless, it is because of both the time-stepping nuance mentioned above, as well as the fact 
that some systems require maximal resolution of highly localized fluctuations in the concentration 
in order to be well-suited to the particular model system, that we choose to model our quiescent 
reactor systems by way of an explicit LDG numerical scheme. We also note that this particular 
numerical scheme has the advantage of being relatively easily generalized to advection dominated 



2 Formulation 



5 



compressible regimes, in which case capturing numerical shock profiles becomes a concern, and is 
often more easily dealt with in the explicit formulations. 

More specifically, we introduce a generalized approach to modeling quiescent reactor systems, 
the theory of which is largely inspired by Ref . [321 123 EU ESI E2] • In §2 we provide a formulation of 
the model problem, then develop the temporal discretization and numerical method for performing 
a predictor multi-corrector over the chemical modes of the system. We proceed by showing a fairly 
standard discontinuous Galerkin spatial discretization, and discuss in some detail the iterative 
methods used along with the temporal mode splitting. 

Let us also note here that a number of very nice numerical approaches to reaction-diffusion 
systems already exist in the literature. In addition to the very nice operator splitting methods in 
the temporal space that employ the Strang method formalism [201 EI] and its entropic structure, 
Petrov-Galerkin (SUPG) residual methods have been applied [31] . fully adaptive finite volume 
(multiresolution) methods have been proposed [57], in addition to compact implicit integration 
factor (cIFF,IFF,cIFF2,ETD) methods over adaptive spatial meshes [32], and particle trajectory 
based methods [IT] , in addition to the stochastic methods dealing with substantially more com- 
plicated molecular scale data [5] |27] . Moreover, exponential convergence results have been shown 
for /ip-adaptive reaction-diffusion systems [421 EH1 EI] where boundary layer data must be retained 
in order to achieve full convergence. In this context we introduce the first — to our knowledge 
- spatially dimension independent /ip-adaptive operator splitting SSP RKDG predictor multi- 
corrector scheme for fully generalized reaction-diffusion systems of equations with functionally 
dependent parameters (e.g. &(a)). 

In §3 we derive the exact entropy relation satisfied by the system, which is borrowed and 
extended from the regularity analysis of [31] , then applying this entropy functional in order to de- 
velop an Zip-adaptive scheme that is fully entropy-consistent — which is to say entropy-preserving 
and bounded — relying only on the global gy^ +1 and local pS^^ ^ entropy densities, as well as 

the local change in the density of the entropic jump J?^<^ e ■ 

Finally, in §4 we present some example applications (that were developed in part using a 
C++ finite element library [0]). We address the complication arising from reactive/diffusion 
dominated/limited regimes explicitly, where we provide four example applications, one which is 
strongly reaction dominated in some areas and diffusion limited in others (a set of fast hypergolic 
combustion reactions), one which is strongly diffusion dominated in some areas and reaction limited 
in others (a set of gas-phase alkyl halide atomic transition metal reactions), and one that oscillates 
between all four regimes (autocatalysis in excitable media across oscillating reactions). Finally we 
utilize an equilibrium system in order to demonstrate the standard and expected error convergence 
results for the method. 



§2 Formulation 

§2.1 Governing equations 

Consider the stationary reaction-diffusion system, which in chemistry and chemical engineering 
contexts generalizes our notion of a quiescent reactor, comprised of % — 1, . . . , n, species in N = 1, 2, 
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or 3 spatial dimensions, satisfying the system of equations: 

dtcii - V x ■ (&i{a)V x ai) - ^J(a) = 0, 

tn j, n 



(2.1) 



with initial-boundary data given by 

Oi(t = 0) = a i)0 , and a^a^ + V^a^ (6* • n + Q • r) - # = on dQ, (2.2) 

taking arbitrary functions a, = dj(t, Xb),bi = bi(t, Xf>), Ci = Ci(t, Xb) and $ = gi(t, Xb) restricted to 
the boundary, where n is the unit outward normal and r the unit tangent vector at the boundary 
dQ. Here on is the concentration of the z-th chemical constituent, the 3>i{pt) are the inter-species 
diffusion coefficients which form an iV x iV x n tensor that characterizes the directional dependence 
on the concentration and its gradient, while v ir G N and v\ r G N are the forward and backward 
stoichiometric coefficients of elementary reaction r G N (if r is not elementary then u( r , v\ T G Q), 
and kf yT , kb, r £ K are the respective forward and backward reactions rates of reaction r. 

More precisely, we are interested in systems comprised of n distinct chemical species OJtj (re- 
calling that each constituent's corresponding concentration is given by a») indexed by r G SH 
elementary chemical reactions, for £H C N where 

4rM 3 === £ <rM k , for r G K, (2.3) 

and where fc/ iT . characterizes the forward rate of reaction r, and kb, r the backward rate of reaction 
r. The indexing sets & T and ^ r are the reactant and product sets & r , 3P r C N for reaction r, 
respectively. The the stoichiometric coefficients Vj r G Q + of the products and reagents, 

and for elementary (or fully reduced form) chemical reactions are positive integers v^ r G Z + when 
i G ^ r or j G since in elementary reactions atoms may only react as absolute entities (which 
is to say in whole number quantities). Furthermore, all chemical constituents of the flow are either 
reactants or products, where inert species may be viewed as the product of a unimolecular reaction 
denoted rj, where z// ■ = v\ rj for all j G 0P ri and k G M VI , such that we may view sri C £? r 
and C .y r . 

Let us proceed by introducing the following n dimensional state vectors: 



a. = (an, . . . , a n ) T , 3> = (0 x (a), . . . , ^(a)) T , *f(a) = K(a), . . . , <(a)) T , 
and additionally defining the "auxiliary variable" cr, such that using stf = stf (a) we may recast 



(2.1) as the coupled system, 

ct t -V x ■ - srf = 0, and cr - V x a = 0, (2.4) 

where we denote the spatial gradient as, V ' x cx = ^f =1 d Xi cx. 



Then we solve (2.4) by employing a predictor multi-corrector method coupled to an RKDG 
scheme which is solved over a reaction mode time-splitting. Let us first describe the solution 
vaguely in terms of three basic steps. In the first step, we use a predictor multi-corrector method 
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to solve (2.4), where we exploit the fact that a partially decoupled version of the reaction source 
term may be solved analytically in order to generate a prediction of the concentrations (aj's) at 
each timestep. This predictor is then corrected by way of a fixed point iteration. In the second 



step, we solve the components of (2.1) in the usual DG formulation, by integrating against test 
functions in space and determining local approximations for each of those terms, respectively, were 
we use an arbitrary order time integrator. Finally, the third step simply requires determining the 
"fast" and "slow" modes with respect to the reactivity of the system, such that for some smallest 
positive n G N we may iterate our solution until nAtf = At s , where we then proceed with the 
same procedure over all the reacting modes (i.e. for each oti G a). The following sections are 
devoted to deriving this methodology. 

§2.2 The predictor multi-corrector 

First notice that the reaction term &/ n = s^ n (a) may be viewed as the source of a proliferating 
set of nontrivial numerical difficulties. That is, not only is it well known that srf 11 may cause 
numerical instabilities due to the varying "fast" and "slow" timescales discussed above, but due to 
the presence of nonlinearities arising from the stoichiometric coefficients z/j ir , the s^ n (a) term is 
responsible for generating a coupled system of n highly nonlinear first order ordinary differential 
equations. Thus, in order to formulate a computationally realistic numerical method for solving 



our system (2.5) over some modest (yet realistic) number of constituents n, we find it necessary 
to employ the following linearization. 

Let us first denote the vector ^(/3j, a) as 



rem \ i=hi¥=i 3=ij¥* 



u b . 



where a is treated as constant for all j ^ i. In other words, we wish to think of 8#\{fii, a) as the 
mass action vector such that all but is treated as temporally inert. 



Then using this notation, we proceed by considering the system of equations (2.4) and dis- 
cretizing in time, such that at time t n+l we are interested in solving the semi-implicit system of 
equations: 



a n+l _ a n 

At 



V x -(® n V x a n ) + ^(a n+ \ot n ), and <r n = V x a n . (2.5) 



The predictor term £/(ct n+1 , ct n ) is predicated on the notion that we may "predict" the ap- 
proximate value of the coupled reaction rate of each species VJHi at time t n+1 by decoupling the 
n-th order system of first order ordinary differential equations containing the nonlinear v{ and 



v\ r factors (as seen in (2.1)) by simply using an analytic rate law derived with respect to the value 
at the previous timestep t n . 

That is, for the i-th molecular constituent 9Jtj we predict its concentration oli at time t n+1 , in 
either the reactant or product well (i.e. in the reaction coordinate representation), by analytically 
solving the following first order ordinary differential equation, 



d t dii = ^(oi^ol). 



(2.6) 
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We solve the integrated rate law in time over a discrete timestep At = t n+l — t n , such that we 
obtain an analytic form for each <5™ +1 in terms of the solution at the previous timestep t n treated 
as constants (again as denoted cx in the term gtf («j, a) for all species of index j ^ i). We refer to 



these solutions (of (2.6)) as the "predicted values" of a™ and denote them either componentwise 



~ n+l 



by a™ or in vector form by a 

Now, of course, the term £/(at n+1 , cx n ) abstractly represents the rate of change of the concen- 
trations «j's at time t n+1 , while the vector ct n+1 will be an averaged form of the predicted total 
concentration at time t n+ . Thus, in order to find the rate of change predictor in the residual 
representation s$ (ct n+1 , ct n ) we simply find the formal difference, 



- n+l _ 

( Xf )• (2.7) 

Next we implement a fixed point iteration corrector over i — 1, . . . , t iterates in order to provide 
convergence in the solution a n+1 . That is, for fixed timestep At we denote the i-th corrector of 



(2.5), implementing our predictor to find a corrector over each iterate i = 0, 



(a n+1 ) i+1 =a n + At (V x ■ (® n o n ) + J((a n+1 ) t+ \ c* n )) , (a n+l )° = (a n+1 )°, 
(6t n+1 ) i+1 = [ +(Q > ) , (T n = V,a«. 



It is important here to recall that the (a n+1 )'s are explicitly determined by the derived analytic 



solutions to (2.6), which depend on the specific reaction system. It is also worth noting that (2.6) 
is solved iteratively in the sense that dta^ 1 = {a n+l Y) is integrated over At to form 

what is the analytic rate law of the new predictor at the (i + l)-st iterate of timestep t n+1 . Note 
here that we have chosen a splitting between the explicit diffusion terms and the semi-implicit 
reaction terms, which we have found (by trial and error) to be the correct splitting to maximize the 
robustness of our method, with respect to both fast and complicated multistable reaction regimes 
(see §4 for examples). 

The endpoint of the iteration i is chosen in tandem with the following bound on the "relative 
change" of the iterated corrector componentwise in j, (i.e. the component (a" 



3 



(a] +1 y-i\\ L ~ (n) 



< C, (2.9) 



for a judicially chosen constant C (e.g. see §4 where C G {10 , 10 }). In fact, we set a slightly 



stronger condition than (2.9), after spatial discretization in the discontinuous Galerkin setting, as 
we provide componentwise convergence in the above sense with respect to each quadrature point. 
Note that we must converge componentwise in each chemical species DJlj, since the relative orders 
of concentrations over !H may substantially vary, such that while the rate limiting products may 
readily converge, the coupled ancillary products may oscillate wildly, and vice versa. In order 
to circumvent this pathology, we demand convergence componentwise globally in space for each 
timestep. 
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§2.3 Spatial discretization 

Now let use discretize in space. Take an open Q C ffi. with boundary dQ = T, given T > such 
that Q T = ((0, T) x fl). Let ^ denote the partition of the closure of the polygonal triangulation 
of Q, which we denote i\, into a finite number of polygonal elements denoted Q e , such that 
= {O ei , Q e2 , . . . , fi e „ e }, for tie £ N the number of elements in Qh- In this work we define the 
mesh diameter h to satisfy h = miny (djj) for the distance function dij = d(xi, Xj) and elementwise 
edge vertices Xi, Xj G dfl e when the mesh is structured and regular. For unstructured meshes we 
mean the average value of h over the mesh. 

Now, let Tij denote the edge shared by two neighboring elements Q ei and £L e . , and for % G I C 
Z + = {1,2, . . .} define the indexing set r{i) = {j G / : Q ej is a neighbor of fl e ,}- Let us denote 
all boundary edges of VL £i contained in dVL h by Sj and letting I B C Z~ = { — 1, —2, . . .} define 
s (0 = {i e : is an edge of such that = Sj for Q e . G £\ when 5j G <9fi e4 , j G Jb- 
Then for H,; = r fi) U s(i), we have 



<9fi ei = [J r y , and <9f^ncK^ = |J Tij. 

We are interested in obtaining an approximate solution to Z7 at time t on the finite dimensional 
space of discontinuous piecewise polynomial functions over Q restricted to given as 

S*(Q h , <%) = {v : v\n H G ^ d (fi e J Vfi 6i G ^} 

for & d (p, ei ) the space of degree < d polynomials over Q £i . 

Choosing a set of degree d polynomial basis functions iVj G & d (Qi) for / = 0, ... ,p we can 
denote the state vector at time t over fi/j, by 



a h (t, x) = J2 "ft*)^ (*)» Va; G ^ 



where the A^'s are the finite element shape functions in the DG setting, and the ckJ's correspond to 
the nodal unknowns. The finite dimensional test functions ip h , $h G W 2 ' 2 (Qh, ^h) are characterized 
by 

d d 

tp h (x)=J2vi N iW and q h {x) = ^\Nt(x) V* G ft, 

Z=0 Z=0 

where and <^ are the nodal values of the test functions in each fl e ., and with the broken Sobolev 
space over the partition & h defined by 

w k > 2 (n h , sr h ) = {v : «,„ g w k > 2 (n e j VQ ei g sr h }. 



We thus multiply (2.8) by the test functions <p h and and then integrate locally over elements 



Q £i in space, defining global scalar products, (ajj, b/Jng = Xln e e% in a ^ : sucn that we 
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obtain the system: 
1 

At 



- oc\ <p h ) = (V, • (^»), <p h ) Qg + K««" ' T ' a-;, ^ 



n+l\i+l ?2 \ 



(K +1 ) ,^) n =((«" +1 ) ,^) 



((«" +1 ) i+1 ,^) fic 



(2.10) 



(^.«fc)n a -(V,a",Cfc)n o = 0- 



We proceed by approximating each term of (2.10) in the usual DG sense. That is, we approx- 



imate the first term on the left in the first equation in (2.10) by. 
1 



ra+l\i+l n 



1 



(2-11) 



Now, let riq be the unit outward normal to dfl ei on T^, and let <£>|r y and ( -P\r ji denote the 
values of ip on considered from the interior and the exterior of Q ei , respectively. Then the 
second term of the first equation in ( |2.10[ ), after an integration by parts, yields, 



(V x -(@W),cp h ) ng = Yl I V x -(cp h @ n cr n )dx-(@W,V x <p h ) ng , (2.12) 



such that we approximate the first term on the right in (2.12) using a generalized viscous flux 
(see Ref. [3]) across the boundary, where upon setting ^ = ^(cr^, tp h ), we approximate 



TL \ TL I Tt I Tl I \ i J 1- 

r fc|r«,0"hl r ji; , Q;^ | r 8j , a h | rj . , riy ) ■ <p h \ T .. 



r N 



(2.13) 



iesf(0 

while the second term in (2.12) is approximated by: 



J? = ^*K, off, <p h ) = {9*1 <pl) » (W, 



(2.14) 



The reaction term in (2.10), which is dealt with using the predictor multi-corrector, is projecting 



into the basis in the obvious way, 



,<Xhh<Ph 



" ) ),<Ph 



(2.15) 



where the second and third equations in (2.10) are approximated componentwise in the usual 
sense, simply setting: 

((" n H +1 ) ,<P h )n «((« B+1 )°^ fc )n fl . 
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Finally, for the fourth equation in (2.10) a numerical flux is chosen which satisfies: 



ire* JTu 



(2.17) 



where ^ / a(a$\ Tij , a£| rji , ?/»|r tf , n -)<E « 

jeS(i) JVi i jeS(i) T » s=l 



§2.4 Formulation of the problem 



Combining (2.11), and (2.13)-(2.17) while setting = JTJg.g^ ^i, we formulate our approximate 
solution to (271) via ( 2.8[ ) which by construction may be stated as: for each n > 0, C E [0, 1] and 
£ > 0, find the pair (ck£, <t^) such that: 



The Predictor Multi-corrector DG Solution 



a h eC\[0,T)',Si), <r h eS, 



/c 



b) (( 



a 



ra+l\i+l 
ft J 



a 



hi 



At(^ + ^) + (^(K +1 ) i+1 ,0,^ 



C ) (K +1 ) ,^) =(k +1 )°,^) ! 



(a n+1 Y 



[a 



<Ph 



a 



j,h ■ 



<Ph\\L°°(n) 



<C7, V(^,j)e(n ei ,^), 



(2.18) 



/) J2f (&, eft * fc , «£) = 0, 

#) 0^(0) = U h ct , V^a^O) = n h a . 

Here, Hh is a projection operator onto the space of discontinuous piecewise polynomials S^, and 
where below we always utilize a standard L 2 -projection, given for a function / e L 2 (f2 e J such that 
our approximate projection / 0(i e L 2 (f2 e J is obtained by solving, L f oh Vhdx = L fo v hdx. 

The gradient projection 11^ merely approximates the initial gradients using numerical difference 
quotients, for example below we frequently employ the approximate fourth order scheme: 



fjx - 2h) - 8/Qe - h) + 8/Qr + fe) - /(x + 2h) 

Y2h 



Now, it follows that (2.18) is a solution for any reaction scheme of arbitrary order satisfying 



the law of mass action (2.1), whether or not the reactions in ^(a) are elementary, and regardless 



of the reaction order {e.g. mixed reaction orders of arbitrary type that change order during the 
course of the reaction, and fractional order reactions, etc.) 



We conclude using a standard time discretization for (2.18), where we employ a family of 



SSP (strong stability preserving) Runge-Kutta schemes as discussed in [581 EH]. That is, we may 
abstractly represent our ODE in (2.18) by = C(oc), such that previously the first order forward 
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Euler time discretization was assumed (equation b in (2.18)). However, to generalize to a 7 stage 
SSP Runge-Kutta method of order p (which we denote SSP(7, pi)), we simply augment the second 



equation in (2.18) in the diffusion terms by: 

,.0')v+2 



h ■ 



b u ) ((cx^) i+2 , ( p h )) =J2(^ k h + At\ jk £ k ,v h ) forj = l, 



,7 



(2.19) 



k=0 



biii) 



a 



n+1 



OL 



hi 



where C k = C(a k l ) can be viewed as the abstract diffusion operator, and the solution at the n-th 
timestep is given as ot% = ct h ^ t=t n and at the n-th plus first timestep by ct'^ +1 = ot h ^ t=t n+i, where 
At = t n+1 — t n . The order p of the method is fully determined by the choice of the coefficients 
\jk and \jk in the Butcher tableau. More clearly, the abstract operator C k does not alter the 
mass action term containing =2/((q;^ +1 ) 1+1 , ck^) in (2.18). That is, as already explained in §2.1 



the predictor multi-corrector scheme provides an approximate solution chosen with respect to 
a distinct temporal integrator which is strongly dependent on the exact solution of a partially 
decoupled system of ODEs; and is thus taken outside the Runge-Kutta loop. 

It remains to identify the "fast" Atf and "slow" At s modes of the system with respect to the 
form of the equations (2.4 often occurring on substantially different time scales, as illustrated in 
[20] . However, the fast modes and the slow modes separate not only with respect to the reaction 
coordinate (as determined by &/), but also separate with respect to the interspecies diffusion as 



determined by Fick's law, which may also be decomposed into fast and slow moving modes. To 
account for this complication we simply decompose the concentration over "fast" <x.f and "slow" 



a s variables, using ex = (a.f,a. s ) 



Such an operator splitting into fast and slow modes has been thoroughly presented in [20] . 
and is known to lead to a fully well-posed system of reaction-diffusion equations, providing the 
existence of an entropic structure and a partial equilibrium manifold. The entropic structure 
discussed in [20J is quite strong however, requiring for example the existence of a C°° monotone 
entropy functional. Clearly in the context of our approximate variational solution (2.18) this 



constraint is too restrictive (esp. with respect to a discontinuous polynomial basis and with an 
eye towards generalizing in order to easily extend our formalism to the reactive multicomponent 
Navier-Stokes regimes, for example). 

Thus in order to stabilize our method we introduce an exact entropic restriction as outlined in 
§3 below, based on the analytic well-posedness results of A. Vasseur, T. Goudon and C. Caputo 
[ISlEl], which depends strongly on an explicit analytic entropy functional We enforce entropy 
consistency on our solution, which provides for the usual monotonicity constraint on the systems 
entropy, but we further expand this constraint and then utilize the entropy consistent scheme as 
a foundation for a dynamic Zip-adaptive strategy as fully derived in §3. 

We further note that the notion of "fast" and "slow" modes here is made to highlight a quali- 
tative choice, where the physics of the system may, of course, be substantially more complicated. 
That is, for simplicity in our derivation, we have assumed that the rate laws split into no more 
than two distinct sets, of "fast" and "slow;" while there may of course be k arbitrary such sets 
representing k grouped rates each of a quantitatively different order of magnitude. While in some 
physical systems it is essential to neglect the chemical kinetics of reactions occurring on substan- 
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tially different timescales (e.g. neutrino production rates in atmospheric chemistry, etc.), in many 
settings (such as in environmental science, for example) it is important to include reactions oc- 
curring in a number of different phases (i.e. ice, water, water vapour, etc.), which can have a 
large array of different timescales for their coupled rates laws. More explicitly, in standard units, 
common chemical reaction rates can differ in a particular setting and choice of units up to some 
twenty orders of magnitude. 

Thus the solution obtained from (2.18) trivially lends itself to these split time discretizations 
(Atf, ■ ■ ■ , At s ), and thus we numerically integrate over the "faster" variables (i.e. the «j G atf) 
in Atf. That is, for some smallest positive neNwe recurse our solution until nAtf = At s (thus 
slightly restricting the permissible choice of Atf), where we proceed with the same procedure over 
all the reacting modes (i.e. Wcti G ex.). In this way we easily acquire the approximate solution over 
the time split modes. 



§3 Entropy enriched hp- adapt ivity 

We are now concerned with an entropy based Zip-adaptive methodology which does not depend on 
the powerful, though often computationally prohibitive, adjoint problem formulation which relies 
upon the calculation of a posteriori error estimates and minimization techniques (e.g. see [16J). 
That is, our method is potentially cheaper computationally, while still maintaining a rigorous 
foundation based on satisfying a priori entropy bounds. 



§3.1 Bounded entropy in quiescent reactors 

We use a formal entropy inequality discovered in [34j in order to develop a stability analysis of 
our approximate solutions. 

Let us derive the entropy of the system J^r over each reaction r G £H in the quiescent reactor. 
First define the mass action term reactionwise, such that we have 



n 

,b 



3=1 3=1 



Now notice that d t oii = «j9 t (lnaj) such that in addition to multiplying (2.1) by lnaj and summing, 
after integration in x we obtain the relation: 



y / ailnaidx — / \^ In aiV x ■ (&i(a)V x a,i)dx — I y J2j r (lno^ + l)dx = 0. 



d 
dt 

Now we use the observation employed in [31], which notes that regardless of the form of the a^s, in 
a system of elementary reactions there always exist reaction dependent constants (b\, . . . ,b n ) G N n 
for hi such that the stoichiometry satisfies the following linear relation: 



$>»{r = Z>*£.. (3.2) 



i=l i=l 



and thus yielding formally upon integration of (2.1) that 4 Ym=i In °i®idx = 0. 
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Now, we rescale the last term by a constant (lnK eq>r = \n(kf^/k b r )) which corresponds (as we 
show below) to the standard Gibb's free energy of the reaction A r G$, such that integrating by 
parts and passing to the weak form we obtain the inequality: 



d 
It 



i=0 

which we rewrite as: 



n „ n „ 

i=0 Jn i=0 Jn 



(3.3) 



n d f n f n 

— / Q!j(lnQ!j + bj)dx + / Qj^(a)V :r Qi ■ V x ajdx + ^j33(o;) < 0. 

i=l i=l i=l 

That is, the reaction term Yli^( a ) = A r G$ corresponds to the isothermal Gibbs free energy 
of the reaction, since 

n n , h f \ 

= In UkI^<A 

i=i i=i ^ ' 

=-ib{ h 1 f ) ^ f in °? ,r ~ u{ ' r + n_i in ^ 



(3.4) 



i=l i=l 



= £ln A^.r + ^/, r n a?" - k bjT J] «?' r ) In Q r (a) 
= £lnQ r (a) + AG e 
with reaction quotient Q r {ot) given by: 

Qr{<*)= (n«^/n^ r )' (3 - 5) 

\i=l ' i=l / 

and where £ = £(£, a;) is simply a rate-scaled prefactor coefficient. Thus, for spontaneous reactions 
at constant temperature $ (the only interesting case for the quiescent reactor systems of the form 
(2.1)), it follows that A r G$ < 0. More precisely notice that we may rewrite (3.4) as: 



n n , h f \ 

£ ®(a) = - £ $i,r ^ ( a^f^ ) 
i=i i=i ^ ' 



= -h, r [ K eq , \{ off.' - J] o^' ) In ( J] of-- / \{ at 

<0, 



(3.6) 



i=l i=l / \ i=l ' i=l 



since it is clear that we have a term of the form (A — B)(XnA — In B) such that the product is 
always positive. 
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As a consequence we obtain the scalar entropy e5^ i0C = =5^ >00 (a) over the reaction space £K. 
That is, given bounded initial reaction state density Po|re« satisfying 

n „ 

Poorest = y~2 a i(^ na i + b i)dx < oo, 
i=i J n 

where a® = (Xi\t=o is the initial condition, then summing over all reactions r G SH we obtain the 
following inequality on the system for any fixed n number of constituents over an unbounded 
domain: 

, n „ n „ t p 

^moo— sup < 2. / a iO- na i + bi)dx + 2. 2_. / / r jD{ot)dxds 
o<t<T I . =1 Jn re?H i=1 Jo Jn 

+ ^2 / a { 1 S# i {a)V x a i ■ V x ondxds > < P \y r m- 
i=1 Jo Jn ) 



where the first term corresponds to the entropy contribution from the density of states, the second 
to the contribution from chemical energy production in the reactor, and the third to the entropy 
contribution due to the random motion of the system. 



§3.2 Consistent entropy and p-enrichment 

The entropy relation derived above may be used to generate a local smoothness estimator on the 
solution of each cell's interior. Moreover, the entropy is a particularly attractive functional due 
to the fact that first, it is globally monotonic and convex (or concave up to the sign convention), 
and also that it is a functional which approximates the local internal energy of the entire solution 
space in a totally coupled sense. In this way, the entropy functional provides for an easy test 



of whether the full approximate solution (2.18) is entropy consistent in a fully coupled sense, and 
then, if it is, where inside of Q the entropy is most variable. In the setting presented in §3.1, 
we have assumed a noncompact domain, and thus J/^ applies globally in the numerical setting 
to periodic boundaries, or those employing fully transmissive (or radiative) boundary conditions 
{e.g. no forcings on the BCs up to the differential order of the numerical solution). 

That is, in order to derive the global discrete total entropy at any particular timestep 

t k+1 , we may simply integrate in time such that for any discrete t e G {0,t h+1 ] we have: 

= sup S.J2 I <4i^<4 + h)dx 
o<t e <t*+i I ~[ Jn g 

+ V" / / aj l S>i{a)V ^ • V x ctidxds (3.8) 
i=1 Jo Jn g 



,_i Jo Jn Q J 



rem i=l ^ u 

where as above a\ = am =t e. 
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Generally, we solve (3.8) when 3>i[a) is any (possibly nontrivial) matrix (e.g. see §4), such 
that in the numerical setting we must compute the following approximation to (3.7) via: 

^£c 1= sup ( J2 I <4{\n<4 + bi)dx 
o<t' ! <i fc + 1 y—[ Jn g y 

+ J2 / Mcu>L} ( ) V X • V x a\dxds (3.9) 

• =1 Jo Jn g \ a i J 

4-k + l 



12^2 / ® S ( a ) dxds < P 0|VreK, 

re <H i=l ^° 



given a small positive constant L G 1R + where 1{ Q4 >l} is the indicator function over the set 
containing a« > L. 

Similarly, to find an approximation of the discrete local in (t,x) entropy we simply 

integrate over an element f2 e . restricted to t +1 , such that we obtain: 



(3.10) 



such that a^ +1 = 0^=^+1, again with L G M + a small positive constant. Then we proceed by 
defining the local in time change in entropy density over int(f2 e J as satisfying: 

wifa = p - -^Uh) > ( 3 - n ) 

where the cell density is taken as p = |fi e J _1 . 

We use equation (3.11) as an approximate measure of the variation in the local internal energy 
with respect to a fixed volume elements (at timestep t k+1 ) interior, int(f2 e J at constant temperature 

More explicitly, we use (3.11) as a local smoothness estimator over the interior of Q ei in order 
to develop a p-enrichment functional <£ k+1 = £ h+1 {& s (fi^ +1 )) which estimates the local internal 
energy of the element as an approximate measure of the smoothness of the solution, such that 
for ^ Pmax (f2 e .) the maximum polynomial order allowed on any Q e ., and ^ s (Q k .) the present 
polynomial order, we define: 

ek+1 = { ^ + W 1 ) if ( I A P^£ " ^KS H I < »s + l) A (S + 1 < Anax), 

1^-1(^+1) if (jApY^ - AgY^J >pj) A( S -l>p min ), 

where the change in the average global entropy density Agy£ +1 at t k+1 is given by 

AgY k+1 = g - S%) ■ 

The global entropy y£ +1 at timestep k + 1 is defined fully in §3.3, and the global density is 
simply taken as g = \Qg\~ 1 . The adjustable parameter /i s = /x(t s ) is a composite of the range of 
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the entropy change at time k + 1 and a weight t s G (0, 1). That is the function may be written 
A*s+i — '•s+i^ over the midpoint of the range S = 8(p, n , S^t t, 1 ) of the change in entropy 
density 



max ApYj^ - min ApJ^ 



fe+i 



Clearly (3.12) has the effect of using the variation in the entropy change locally to weight some 
fraction of the cells for p-enrichment and the rest for p-coarsening, depending only on how far 
their relative local change in entropy lies from the median of the range. 

In contrast to alternative choices for p-enrichment, this scheme provides for a cogent physical 
interpretation that serves as buttress for the enrichment strategy. Namely, we see that in areas in 
which the relative disorder (i.e. the relative entropy) of a cell exceeds a specified allowed variation 
within the cell itself, then we coarsen our solution, thus avoiding CFL instabilities, etc. Likewise 
in areas of relative order (or stable smoother regions) we readily enrich our solution. 

We have also tested enrichment strategies based on slightly more abstract principles. For 
example, one may simply choose a fraction of elements with respect to the magnitude of their 
relative change in entropy density, or with respect to | ApS^t^ — A@yJ^ | by cell. Likewise, when 
using a hierarchical basis one may use the scheme described in [50] to measure the perturbative 
variation in the higher terms with respect to the L 9 -norm. Many of these alternative strategies 
can lead to stable schemes that effectively "sense" relative energy fluctuations with respect the 



At. However, it should be noted as a word of caution that (3.12) is particularly well-suited for 
naturally avoiding the observed phenomenon of bunching in the local variational space. This 
bunching of the solution often leads to a flickering of enrichment /coarsening of a substantial 
number of elements taking values close to the "center" of the chosen discriminating parameter (e.g. 



|ApJ^Q* — Agy£\l e | in (3.12)). This behavior is a potential source of debilitating inefficiency 
in the scheme, and can be difficult to isolate without recourse to an entropy formalism, or similar. 



§3.3 The entropic jump and fyp-adaptivity 

As already discussed, the global entropy formulation from §3.1 is predicated on noncompactness 
of the space Q, while in general we are interested in more complicated boundary formulations; 



in particular any boundary condition satisfying (2.2). In this more general setting we see that 



equation (3.7) by way of the divergence theorem becomes: 



{n „ n „ t p 

/ / ai(\n ai + bi)dx + y^y^ / / D(a)dxds 

J J a^ 1 3ii(a)V x OLi ■ V x ctidxds | (3.13) 

< Po\wem + ^2 / ( m + bi)^i(a)V x ai ■ ndSds. 
Jo Jon 



E 



8=1 
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Thus, as before, the discrete approximation to (3.13) simply yields: 



= sup af(lna| + 6j)dx 

o<t<<t*+i \~[Ja g ' , 

+ / / M^>l} ( ) • V x afdxds 



i=l " u 



(3.14) 



+ ^2^2 / < P |Vre« 

reSR i=l ^° ^ 

n r tk+1 r 

+ J2 / l{a i >L}(lna? + 6 i )^i(a)V a! aJ -ndSds. 

i=1 ^0 ./d^g 

Further, notice that for the appropriate choice of boundary conditions both = J^oo and 



Now, in the local approximation it is clear enough how to reformulate (3.14) over cells such 
that we obtain a local approximation to the entropy in the neighborhood of the cell. However, 
for the case of /i-adaptivity we are more directly concerned with the local jump in entropy across 
the neighboring cells, since it is these jumps which serve as a proper diagnostic probe for stable 
/ip-adaptivity (e.g. see (TTJ HH HQ]). Thus we define the local entropic jump at time 

t k+1 by 



tk + 1 



fi$L. = E / / Mol^o* + h)%(a)V x al ■ ndSds, (3.15) 
■ i=l Jo Jan.. 

such that the density of the change in the entropic jump pA ^J^ e is given to satisfy 

P^St&t = P (Xn. " A) • ( 3 - 16 ) 

We proceed by estimating the approximate flux of the internal energy of the system by con- 
structing the Tt-adaptivity functional 21 = 2l(e5^/(fig +1 )) where the mesh triangulation 3?h at 
time t k given by h = h(t k ,x) is refined to level h' = h(t k+1 ,x) — that is, we isotropically re- 
fine to h/2 in each spatial dimension — over cell Vt ei at time t k+1 . Similarly we may unrefine 
^fc+i = 2l(«^o(n* +1 )) to level h = h(t k+1 ,x) — that is, we isotropically coarsen to 2h in each 
spatial dimension. 

For example, in dimension N = 2 the refinement would take a quadrilateral parent cell Q ei 
and split it into four child cells Cj, while a coarsening would take four child cells denoted Cj and 
merge them into a single parent element Q ei . Thus depending on the evaluation of 21, we obtain 
the full /i-adaptivity functional: 

o,/m / M^ 1 ) if (ipA^Jt - gAj? k+1 \ > rj v ) A (s + 1 < h mgx ), 

21 : — < / v _ x ' (3.17) 

\& h Mt X ) if - ^ A X +1 | < Vh ) A (S - 1 > ^nin) VC„ 

where /i max , ft-m in correspond to the maximum and minimum refinement levels, respectively. Here 
again, the density of the global change in the entropic jump is given such that: 
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where q is the same as in §3.2. Also as in §3.2, the adjustable parameter rjh = vi v h) is again 
defined over the range of the change in the entropic jump ip = ip(p, J?Jl{ ne , ), and is given 

by Vh = Vhip such that 

V> = maxpAjr*+t ei - rninp^Jv 

and Vh € (0, 1). 

It is further interesting to note that the local change in the entropic jump A J?^^ is indepen- 
dent of the reaction entropy at time level t k , and ends up depending only on the reaction coupling 
from the earlier timesteps as well as on the present states diffusivity. That being said, it is clear 



that just as (3.12) in §3.2 effectively p-enriches the solution based on the local physics of the 



system, here, we find that (3.17) has the effect of flagging elements with a high relative change in 
their entropic jumps for h- refinement, and those with low relative change in their entropic jumps 
for /i-coarsening. That is, in areas where the entropy is changing dramatically across the elements 
boundary, we refine. However, when coarsening, we are presented with the additional constraint 
denoted: VCj. That is, by VC^ we simply mean that in order to actually coarsen a parent element 
fl ei comprised of j children elements DjCj = fl ei , each child Cj must be independently flagged 
for coarsening. In other words, all children of an isotropically refined element Q ei must contain a 
coarsen flag at time level k + 1 in order for the parent cell to ultimately be refined at time level 
k + 1. For more details on this isotropic refinement strategy we direct the reader to [7]. 

Finally, we couple the /i-adaptivity functional %l k+l to the p-enrichment functional <£ k+1 such 
that /i-adaptivity is always preferentially chosen over p-enrichment. That is, clearly the cell lo- 
calized entropy and its corresponding entropic jump are strongly coupled by virtue 



of (2.1), but in order to avoid numerical instabilities caused by erroneously p-enriching relatively 
inert cells experiencing high entropic fluxes entering through its neighbors, we evaluate the simple 
kinetic switch functional = ^ +1 (2tf +1 , <£ k+l ) determined by evaluating: 



A' 



fe+i 



' a + 1 if & h , (n^ 1 ) a &> s a ,y k+1 , 
oi k+1 if ^h'i^ef 1 ) a & s+ \ri k + l ) a y k+1 



A £fe +1 . f ^(ftfc+i) A £aw(Qj+i) a y k+ \ (3.18) 

ok-\ 



2tf +1 A <B k+1 if ^hol^et 1 ) A ^ s+ \^ k et l ) A 



otherwise, 



whereby we are able to stabilize these spurious quiescent instabilities, and yet still maintain the 
entropy consistency of the scheme. 



§4 Example Applications 

We address several example applications below, and note that all examples in this sections were 
given reflecting wall boundary conditions, which is just to say that scalar boundary values are 
determined by their values on the interior at the boundary ajO;j&|an e . — <*i\m<..i and the gradients 

1 3 

are reflected with respect to the normal direction bi(T i)b \gQ e = — (Ti\on e ., where c, = 0. 
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§4.1 Reaction dominated hypergolic kinetics in ID 

We choose as a one dimensional example the nuanced problem of the reaction dominated — as 
well as diffusion limited - - regime in hypergolic kinetics, corresponding to when the rates of 
the reactions r in !H occur on substantially smaller timescales than the corresponding diffusivity 
@i of the system. In the proper context (such as in a stable subdomain Q d= Cl C f2 of a 
homogenized combustion chamber [23JH7], or in the propagation of deflagration flame fronts [T2] ) 
such a combustion system may be approximated by the reaction-diffusion equations, and to a low 
order approximation — where the compressibility of the fluid may be neglected — the quiescent 
reactor regime may be utilized to model the resulting gas phase dynamics of the chamber. 
Along these lines let us consider the second order hypergolic ignition reaction, given by 

/ i / k l 6,6 



comprised of a reaction between monomethylhydrazine a\ = CH 3 (NH)NH 2 and nitrogen dioxide 
«3 = NO2, to yield the exhaust radical a<i = CH 3 NNH 2 and nitrous acid «4 = HNO2 such that 
the stoichiometry satisfies: 

CH 3 (NH)NH 2 + N0 2 — ^ CH3NNH2 + HN0 2 , (4.1) 



with a forward reaction rate (see Ref. [H]) of kf = 2.2 x 10 11 e _5900 / Rl? cm 3 (mol-s) _1 , where d is the 
constant temperature subdomain Q in which the reaction occurs and R is the ideal gas constant. 
Being far from equilibrium, we can readily neglect the back reaction, such that for a 1 ,a 2 ,«3 



and a 4 we see that (2.6) leads to the coupled system: 

dtoii = —kfcticts, <9t«2 = kfaiac3, dta^ = —kfa^ai, dtot^ = kfaia3, (4-2) 
which upon integration over a discrete timestep At = t n+1 — t n , taking a™ +1 = ai(t n+1 ), yields: 



(4.3) 



In the case of combustion reactions a mass transfer correction factor h m G K is often included, 
and is in fact necessary in order to stave off the effects of catalytic volume expansions, where 
again the implicit assumption is that there exists a local subdomain Qq of relatively homogeneous 
reactivity, over which the rate constant kf may be effectively averaged (e.g. see Ref. [72]). This 
correction-based formulation leads to a stable form of the mass action, such that simply replacing 



kf in (4.3) by the addend (kf + h m ) provides the vector cx n+1 from which we may easily compute 
J(a n+1 ,a n ) from Q. 

It remains to identify the mass diffusion coefficients ^(a, 1?). The most straightforward way 
of choosing @i(a, 1?) is by simply setting them equal to empirically determined constants, in which 
case one may make the additional assumption of pure diffusion, such that interspecies diffusion 
3tij(a, 1?) — for a counter example see §4.3 — is neglected. In this case the ^ reduces to a diagonal 
matrix with positive constant entries, ^ G M + , and setting the transmissive boundary conditions 



Ul\ Kjl = Efft|jCij, we may proceed to solve (2.1) for (4.1). 
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Index 


Species 


Mass Diffusion/cm 2 • s 1 


a,\ 


CH3N2H3 


~ 5.9 x 10~ 6 (1) 


a 2 


CH2N2H2 


~ 6.0 x lO" 6 


a 3 


N0 2 


~ 9.0 x 10~ 2 (2) 


«4 


HN0 2 


~ 1.2 x 10" 2 (3) 



Table 1: All approximate values taken at STP since at constant pressure ^ scales sublinearly 
with •& (e.g. see §4.3). (1) was measured via chronoamperometry as shown in Ref. [52J, (2) was 
determined via the single component Chapman-Enskog experimental fits in Ref. [44J, and (3) 
was calculated using diffusion denuders in Ref. [TO]. The remaining coefficient was adapted using 
relative magnitude arguments from simple collisional theory [36J (viz. ^ oc a~ 2 the molecular 
cross sectional radius and §4.3). 



The results are presented in Figure (|TJ using the diffusion constants from Table 4.1, where the 
initial conditions are given by: 



011,0 



= 0.49 exp 




= 0.49 exp 



(x - 40) 
800 



«2,0 



= 0:4,0 = 1 x 10 



-5 



We use a Bassi-Rebay form of the LDG based viscous flux (see [U EH]), and a mass transfer 
correction factor of h m = —2.19 x 10 10 e~ 5900 / Rl? cm 3 (mol • s) _1 in order to rescale kf to an effective 
value of (kf + h m ) = 1 x 10 4 e~ 5900/ ' Rl? cm 3 (mol • s) _1 . Then our slow (diffusion) modes are chosen 
such that At s = 1\<ls and our fast (reaction) modes to satisfy Atf — eAt s for e G {^j, l}. The 
convergence bound from (2.9) is taken as C = 1 x 10~ 8 over both fast and slow modes, and the 
average number of convergence steps is ~ 25 where i = 50. 

As is clear from figure [TJ running the reaction with the fast modes has a substantial impact 
on the solution, but only up to a point. That is, once a certain timestep is reached with full 
convergence, little is gained (in a relative sense) by iterating with respect to additional temporal 
refinements. Here, as is clear, the diffusion is negligible on the timescale of the simulation, and 
has nearly no effect on the solution over 2 [is. 

We also note that while the addition of both the mass transfer correction h m and the fast 
reaction time modes Atf makes solving such problems viable in the reaction-diffusion setting of 
(2.1 ), most such problems — certainly most highly reactive gas phase or combustion reactions - 
require the addition of advective fluxes in order to accurately approximate the relevant dynamics 
of the system, which are determined by complicated couplings between the mass, momentum 
and energy conservation equations (e.g. the Euler and Navier-Stokes equations), not to mention 
to substantial importance that turbulence plays in these types of combustion regimes. We shall 
discuss these systems in more detail in the sequel to this paper. 
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Figure 1: The top graph shows the solution for h = 0.5, At s = and Atf = At s /50. The 
bottom left shows the absolute difference map between the a's solved using Atf = At s /50 and 
those solved using Atf = At s , while the bottom right shows the absolute difference map between 
the a's of At/ = At s /50 and At f = AtJlO. 

§4.2 Diffusion dominated alkyl halide gas mixture 

Consider the following diffusion dominated and reaction limited gas phase reaction of ffuoromethane 
CH3F, with free gold cation Au + as discussed in [HJ EH] and whose primary reaction channel fol- 
lows the termolecular addition reaction in the presence of He gas: 

CH 3 F + Au + . kfl • Au + CH 3 F (4.4) 

where the helium is an important third body in the gas phase reaction pathway, and is further 
coupled to the bimolecular elimination reaction (neglecting the hydrogen halide formation): 

CH 3 F + Au + . kf2 • AuCH+ (4.5) 

The forward reaction rates are taken from [76], and are given to satisfy kf l rj 8.9xl0 _12 cm 3 /mol- 
s, and kf 2 rj 3.4 x 10 _12 cm 3 /mol • s. The backward rates are simply derived from the equilibrium 
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Figure 2: On the top we show the N = 2 solution using five levels of h refinement given p fixed at 
linears, where the contour is coloured by max x o; on the left and by the S from §3.2 on the right, 
with v h = 0.04 and h G {1/16,1/32,1/64,1/128,1/256} . At the bottom we have 3 levels of p 
at fixed h = 1/64 on the left, and on the right the corresponding values of ip from §3.3, where 
l s = 0.04. Each graph is shown at T = 1.3 minutes. 



constant K eq approximation, using the standard isothermal Gibb's free energy change for the 
reactions AGf = —Rd\n.K eq) as measured in [76] are given as K eq = 2.63 near standard state 
(p = 1 atm, $ = 295 ±2K), such that we can simply assume that = k^/Keg and k^ 2 = kf 2 /K eq . 

Then adopting the standard — though here prolix — notation for a\ = [CH 3 F], a 2 = [Au + ], 
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a 3 = [Au + CH 3 F], and a 4 = [AuCHj], we obtain: 

d t [CE 3 F] = A; fel [Au + CH 3 F] + fc 62 [AuCH+] - (k fl + %)[CH 3 F][Au H 
d t [Au + ] = A; fel [Au + CH 3 F] + fc b2 [AuCH+] - (k fl + fc /2 )[CH 3 F][Au + 
d t [Au+CH 3 F] = (k fl + fc /2 )[CH 3 F][Au+] - k bl [Au+CH 3 F] , 
«9 t [AuCH+] = (% + %)[CH 3 F][Au+] - fc 62 [AuCH+], 



(4.6) 



where [He] is the inert reactive bath. Now, as discussed in §2, it is not difficult to explicitly solve 
this system of first order ordinary differential equations in terms of the previous timestep. That 
is, for an ODE in w of the form w' = C\w + Ci (as is each of our constituents in 4.6) we have the 
general solution over At, 



(4.7) 



Thus letting f(t n ) = f n for our coupled system (4.6), and solving each first order ordinary 
differential equation, we arrive with 



[CH 3 F] 



71+1 



k bl [Au + CH 3 F] n + k b2 [AuCH+] n 



(k h + k h )[ku + Y 



+ e -(*/ 1 +*/ 2 )[Au+rA t / [CHsF] 



fc 6l [Au + CH 3 F]" + fc fc2 [AuCH^ 



(k fl + k h )[A^ 



[Au 



+in+l 



k bl [Au + CH 3 F] n + k b2 [AuCH^ 



(k h + k h )[CR 3 F}" 

+ e -(^ 1 +*/ 2 )[CH 3 F]»At |jA U +] n 

[Au + CH 3 F]" +1 = ±- ((% + %)[CH 3 F]"[Au + ]") 



fc 6l [Au + CH 3 F] n + fc 62 [AuCH^] n 



(4.i 



+ e~ k »i At ( [Au + CH 3 F] n - ((% + %)[CH 3 F] n [Au + ]") 



[AuCH+]" +1 



-^((% + A; /2 )[CH 3 F]"[Au + ]") 



<b 2 



_|_ e -fc b2 At 



[AuCH 



1 



2 r- — ((% + %)[CH 3 Fr[Au + r) 



Now these compri se th e analytic predictor solutions from (2.6), and we are without difficulty 
able to construct (2.7) in the mass action operator 

Note that our initial conditions are taken to satisfy Gaussian distributions in N = 2 dimensional 
space: 

[CH 3 F] Q = 6e -((-l-6) 2 +(y-1.6) 2 )/1.5 ; [ Au+ ] o = ^-((,-0.6)^-3.6)^ 

[Au+CH 3 F] = 4e -((-o.4) 2 + (y-o.4)2)/ 2) [ AuC H+] = 6e -((-+i-6) 2 +( 2 /+i.6) 2 )/i.75 ) 

where for simplicity here, we take as a first order approximation that the diffusivity is constant 
with respect to the helium bath f ~ 5 x m 2 -s _1 (for example see the CRC handbook of 
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h vs. p entropy behavior 

474 1 ' ' ' ' i ' p= 




472' 1 1 1 1 1 1 1 1 ' 

10 20 30 40 50 60 70 80 90 100 

Timesteps 

Figure 3: Here we show the entropy consistency of each scheme from Figure [2j where the uncoupled 
p-enrichment scheme shows a slightly offset decay in S?^ 1 than that of the uncoupled /i-adaptive 
scheme. 



chemistry and physics online edition, and compare helium methane/sulfur hexafluoride mixtures, 
etc.). We show some numerical results in Figure [2j where At = 10 s, the initial h level locally 
is h = 1/16 and the highest level of refinement corresponds to h = 1/256 locally. Similarly the 
initial p level is p = 1 locally, and the highest level is p = 3. The multi-corrector is set with no 
fast timesteps here, and with the convergence bound of C = 10~ 12 taking i = 10, where in most 
steps «i is the slowest to converge, and usually does not achieve the C bound before reaching £ 
— often achieving a C-tolerance corresponding to ~ 10 -9 . 

We also note that both solutions shown in Figure [2] — the uncoupled /i-adaptive and the 
uncoupled p-enrichment schemes — are entropy consistent as displayed in Figure [3] up to T ~ 17 
minutes using timesteps of At = 10 s. However, this example is quite simple, and we have not 
shown the full /ip-adaptive results. Let us now address a more complicated system: the Belousov- 
Zhabotinskii reaction. 



§4.3 Autocatalysis: the BZ reaction with S?ij(a) 

Here we consider a mixed regime, where the reactivity and diffusivity demonstrate a complicated 
and subtle interplay between the reaction and diffusion limited and dominated regimes, respec- 
tively. That is, the classical reaction-diffusion known as the Belousov-Zhabotinskii (BZ) reaction 
is known to demonstrate reactive oscillations in time, and comprises what is known as an excitable 
medium; which is simply a medium whose propagation is nonlinearly constrained by a dispersion 
limited (visible to the eye) refractory period caused by local gradients in the rate limiting reagent. 
In this sense, the non-equilibrium thermodynamics of the BZ reaction illustrates a nice example of 
a reaction regime that oscillates between a diffusion dominated process and a reaction dominated 
process. 

The generalized chemical kinetics of the system are characterized by the following coupled set 
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of chemical reactions: 



BrOg + Br" 
HBr0 2 + Br" 
BrOg + HBr0 2 
2HBr0 2 
B + M 



«/2 



/ ! 



t: /,l 



S /5 



HBr0 2 + P 



2P 



2HBr0 2 + 2M 



BrOg + P 



cBr 



(4.9a) 
(4.9b) 
(4.9c) 
(4.9d) 
(4.9e) 



such that oil = [HBr0 2 ], a 2 = [Br ], a 3 - 
in addition we take a positive constant c G 



c = c(ag) where a 6 \ 



t=o 



[M], a 4 = [BrOg], ct 5 = [P] and a 6 = [B] where 
l + depending on the initial concentration of B (i.e. 

1 l+vg 
4' 2 



and contained in the interval c 6 



whenever dynamic 



oscillations are present in the reaction space 9\. That is, c is considered here as an adjustable 
stoichiometric coefficient, up to an adjustable rate constant kf 5 , which depends on the choice of 
catalyst M and the oxidized (generally organic) species B, where P and BrOg are the generalized 



reaction products. As is the standard convention, we neglect the hydrogen cations in (4.9a)-(4.9e). 



The classical approach to modeling the BZ reactive system is to employ the Oregonator model 
[22] 14*6] . which is a simplification of the FKN (Fields, Koros and Noyes) chemical mechanism 
[29] [30] that relies strongly on the underlying assumption that the products of the reaction are 
decoupled from the dynamics of the oscillation mechanism [67J. It should also be noted that 
alternative model systems do exist, such as the radicalator and GTF models [67J, which generally 
attempt to formulate systems reliant on many more degrees of freedom in a., and to eliminate, 
for example, the adjustable parameter c. These extended models introduce their own problems 
however (e.g. n = 26 and r max = 80), and it remains unclear, at present, whether they model the 
experimental behavior better than the classical Oregonator. 

Below, we expand on the classical Oregonator model to couple all six unknowns of the system 
(including the product well constituent [P]) and derive a new formalism for treating this system 
which utilizes the methods implicit in the scheme presented in §2 and §3. However, expanding 
the Oregonator to dynamically include the bath components leads to a number of additional 
complications in the solution space. For example, it is well known that the Oregonator model 
admits bistable nonequilibrium solutions [22], while here, this bistability is observed, and moreover 
- as we predict from studying the surface of the reaction coordinate (e.g. see Figure [ij — the 
solution in this extended treatment may develop local multistability behavior due to the highly 
irregular reactive surface. However, we also observe a dampening effect noticeable near the limit 
cycle, which is not usually present in the standard Oregonator model, though very nicely agrees 
with the experimental observation in oscillating reaction experiments. This is not surprising since 
local fluctuations in the bath concentration can have large effects on the local magnitude of the 
functional parameters (e.g. @ij(a) and K eq (a)). These observations require closer analysis, and 
generally lie beyond the scope of this paper. Here, we restrict ourselves primarily to introducing 
the model. 

Starting with the quiescent reactor scheme from §2 we partially decouple our kinetic equations 
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for t G (0,T) from (2.1) by solving the following fully explicit integrated rate law: 

d t a,\ =kf 1 a / ±oi2 + kf 3 a^a\ — kf 2 a\Ot2 — 2kf 4 a\ 
d t at2 = ckf 5 ae,oi3 — kf 1 a^a2 — kf 2 a\a2 
dtocs = 2kf 3 aia4 — kf h a%a^ 
dton = kf A a\ — kf 1 a 2 a^ — kj z a\a^ 
d t a^ = kf A ot\ + 2kf 2 aiCt2 + kf 1 a^a2 
d t a e = -k f5 a 3 a 6 . 



(4.10a) 
(4.10b) 
(4.10c) 
(4.10d) 
(4.10e) 
(4.10f) 



We also note that we explicitly solve here for the reaction products B,P and Br0 3 , which are 
conventionally neglected due to the fact that they are each present in excess throughout the 
medium and are known to have negligible effect on the dynamics. However, in our system they 
play a rather central role, as we will see below. 

Now, it is immediately clear that all of the corresponding differential equations (4.10b)-(4.10f ) 
are linear with respect to time, with the exception being (4.10a). That is, the linear equations, as 
in §4, may be solved using (4.7), which yields in the a notation that: 



a 



n+1 



a 



ck h alal 

k h a 2 + kf2 a i 

2kf 3 Q>ia2 

k h< 



+ e 



<k fl a2+k f2 a?)At ( ft n 



n+1 



"3 - 



kfoO-l + kf 2 a™ 
2kf 3 a™a2 



(4.11) 



a 



n+1 



(a 



n\2 



k h a 2 + k f3 a l 



n\2 



k h a 2 + k f3 a l 



a 



n+1 



At(kf 4 al + 2kf 2 (X\(X2 + k^a^x-i) + <y§, and a% +1 = «6 e " 



k f5 a%At 



The first equation, on the other hand, requires a solution to the classical nonlinear Riccati 
equation, which means solving for 

C\uj 2 + C 2 n7 + C 3 , 



(4.12) 



where here we may treat as constants C\,C2,C% G R. Generally we may solve (4.12) by way of 
the fundamental theorem of calculus, such that we first obtain 



1 



ds — t = 0, 



for C4 G K. Then the general solution is simply determined by the sign of the discriminant 
C = Cf — AC1C3 of the polynomial in s from the denominator, such that upon integrating we 
arrive with the general solution: 



zu 



-l/dt - for Cl = 4dC 3 

ICf 1 (tan(i^ 1/2 )e 1/2 - C 2 ) for Cf < 4C 1 C 3 
iCr 1 (tanh(^C 1/2 )C 1/2 + C 2 ) for Cf > 4dC 3 



(4.13) 



Now, in order to find the well-posed discretized version of (4.12) over At™, we note that the 
differential equation in ot\ similarly yields, 



(k fl a™a n 



t 4 ct 2 



+ (fc/3<*4 



k h a^)s - 2k h s 2 ) 1 ds - At n = 0, 
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while the general solution to the differential equation yields 



-2^( ta " h (^ C,/2 ) Cl/2 + C! ) (414) 



where C is a constant that depends on a(0) and ( = Cf — AC1C3 > 0. Now we will be concerned 
with determining given a™ f° r t n+1 = t n + At. Then setting ( n = {C% ) 2 - ACiC^ > 0, and 
letting C\ = — 2kf 4 , C'£ = (kf 3 a2 — kf 2 a%) and C3 = k^oQoQ, we have 



^( ta , h (-f^c-)c- + c ? ) 



but this involves the unknown C and not a™. Then by rewriting the argument of the tanh into 
two groups: 

n+l 



a 



1 ~ 2C ' 

and expanding using 



- (tanh (^±% 2 + fey 2 ) c /2 + c?) 

tanh(x) + tanh(y) \ 



tanh(x + y) 



we see that 



+ tanh(x) tanh(y) / 
tanh (^Cn /2 ) + tanh (f ('J 2 



^1 + tanh (^Cn /2 ) tanh (fC^ /2 ) 



Cn /2 + C 2 n . (4.15) 



Then recognizing that the term from (4.14), 

tanh ( t -^C 1 J^j rearranges to, - (2C ia J + C^Q 1 ' 2 = tanh (^±^ /2 ) , 

such that substituting back into (4.15) eliminates the term with the unknown C in terms of the 
known a™, whereby we arrive with 

1 /-(2d«^ + C 2 )Cn 1/2 + tanh ffCn /2 ) 

a n+l = l _ V 2 / A 1/2 + q 

2C i \ 1 - (2Cia? + C 2 )Cn 1/2 tanh (faj /2 ) " 
so that letting i n = ^JC\G^ and noticing that G\ < we recover the full solution: 

■ (A*(t n a? + C 3 n ) - Cia?) / (At(t n + Ci«?) - Ci) for Cn = 0, 

x ^^^-^w^^)^^^ forCn>0 (4.16) 



2C i ^l-(2CiaJ+C2)Cn 1/2 tanh(^ai /2 ) y 

It remains to address our variable in a diffusion tensor. That is, in a slightly more general 
setting than in §4.1 and §4.2, it is the case that the diffusion tensor is known to obey functional 
dependencies such that, for example ^ = @ij(a, 1?) for $ the temperature at which the reaction 
occurs (thus constant in the isothermal approximation). These dependencies may be determined 
in a number of different ways, namely: (a) they may be determined empirically, (b) they may be 
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determined using basic collisional theory arguments, (c) they may be determined using the Stokes- 
Einstein relation, (d) they may be determined by applying the fluctuation dissipation theory, and 
so forth. In any case, the coefficients ^ may be simply taken as the row sum over the matrices 
%, such that % = YTj=i %• 

We present a simple functional form for the mass diffusion (neglecting concentration gradient 
dependencies {e.g. = %j(V x a,a,$) in [SJ) as discussed and derived in Ref. [151 E31 ES] which 
arises as a natural result of the Chapman-Enskog theory. That is, recall that the j-th constituent 
Wlj is written as a concentration with respect to the specific volume Fixing our <x,-'s as the 
molar fraction of species VJtj, where Mj corresponds to its specific molecular mass, we may write 
the standard multicomponent diffusion tensor as satisfying: 



K ji - K % 




where the cofactor matrices are given by: 



(4.17) 



K 3% 



-1) 



i+j 










■ K 1>n 


















K n ,i . 


Kn,i-1 








(4.18) 



with entries defined by, 



K. 



a; 



+ 



Ma 



E 



if i j, and zero when i = j, 



such that the binary mixtures are set componentwise by way of the reduced molecular mass via, 

[%\ = % j p- l ^/¥iM~+M~)/2MW j . 



The Chapman-Enskog prefactors ^ = %j(ai, aj) (see Ref. [151 EH] for derivation and details) are 
defined in terms of the reduced temperature of the mixture = i?*j(a», a 3 -), a unitless function 

of the first order deviation from the idealized rigid sphere model denoted fl^ , and the cross 
sectional radius <Jij = aij(9Jli,9Jlj) in A such that in the first approximation, ^ = 2.628 x 



g • mol 1 and the 



(viz. 

%A are m cm 



equation 7.4-4 in 

2 „-l 



where is in K, p is in atm, the are in 



For simplicity, we take the deviation parameter to unity to restrict to the rigid sphere 



approximation. Then setting <7j 



Kbl >max , where b ijjmax 



(rj + Tj) such that the T{ correspond 



to the approximate maximum radius of a single molecule of iVli in the rigid sphere approximation 
to the reactive cross section as shown in [36J (where molecular geometries and bond lengths are 
approximated using the ghemical/Mopac suite), and approximating the critical relations in the 
binary mixtures via the weighted sum = J2ke{ij} a k$li as seen m PQ f° r example, where the 
approximate pure critical temperatures are given using the technique developed in 
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Figure 4: We plot (4.16) near a limit cycle when („ > on the left, and the on the right is the 



entropy consistency of the solution from Figure after T = 25 seconds. 



Species 




a 2 


a 3 


«4 


«5 


«6 




.343 


.364 


.428 


.345 


.341 


.428 


n/k 


3.35 


1.2 


6.4 


2.8 


2.3 


6.4 



Table 2: We show the results of the crude approximate values for the reduced temperature and 
the effective maximum cross sectional radius of each pure species 9Jtj. 



Further note that we may rewrite (4.17) using the classical adjoint matrix such that 



(4.19) 



where (K represents the r/'-th entry of the full rank inverse matrix K 1 . 

We take as our example case the system studied in [22], modified to include the functional 
dependencies of f^-. That is we use the ferroin system such that P = HOBr, M = Fe(Phen)g + 
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Figure 5: On the left is the /i-adaptation at T = 3 seconds from initial as a wireframe hexahedral 
mesh, where h G {1/4, 1/8, 1/16, 1/32}. On the right we show the p-enrichment at T = 3 seconds 
with p G {1,2, 3}, where p — 1 is red, p = 2 is green, and p = 3 is blue. The domain is Q = [—5, 5] 3 



and B = M. For full consistency we recast (4.9a)-(4.9e) in the form: 

BrO 



Br" 
HBr0 2 + Br" 
Br0 3 + HBr0 2 
2HBr0 2 
2M 



*/2 



«/3 



t: /,l 



C /5 



HBr0 2 + P 


(4.20a) 


2P 


(4.20b) 


2HBr0 2 + 2M 


(4.20c) 


Br0 3 + P 


(4.20d) 


cBr~ 


(4.20e) 



which means that (4.11) must be augmented in the second, third and sixth constituents, such 
that: 

d t ct2 = ckf 5 al — kf 1 a^a2 — /c/ 2 ai« 2 
<9i«3 = 2kf 3 aia4 — 2kf 5 a\. 
The sixth constituent vanishes here while the second yields: 



(4.21) 



kfj_a2 + kf 2 aii 



ckfja?) 2 



kf 1 a2 + ^/ 2 tt i 



(4.22) 



The third requires solving another Ricatti equation (albeit a simpler one), such that we consider 
the equation a' 3 = C3 — Cia^, where similar to before, we begin with the general solution 



«3 



-^tanhf^ + ^v^a 



(4.23) 
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k h 


k h 


k fi 


k h 


2.5 [H+] 2 M-V 1 


3 x 10 6 [H+] M- 2 s- x 


40 [H+] M'V 1 


3 x 10 3 M'V 1 


0.1 s- 1 



Table 3: The reaction rates used for the BZ reaction. 



and perform the expansion for a n+l 
we acquire 



a(t n + At), such that setting C% = 2kf 3 a™o>2 and C\ = 2k 



a 



n+l 



y/CfC~ 3 tanh ((t n + C) y/CfC^) + tanh (At^/CfC~ 3 ) 
CI 1 + tanh ( (t n + C) y/Cf(h) tanh (At" y/CjC^) 



Then rearranging (4.23) such that 



:« 



we again may eliminate the C, obtaining 



tanh ((t n + C)y^CYa 



a 



n+l 



£L=a« + tanh (At^CfG,) 



£tanh (Aty/CfQs) 



V 



and thus ultimately recovering the full solution: 



' a% + ^gP tanh (At^CfU~ 3 ) ^ 
tanh (Aty/CfC^) 



1 + 



n+l 



^/(2/c /5 ^At + l) 

/c n C'i 

V c n tanh(At x /C 1 "C 3 ) 




Jtanh(At v /C* 1 n C 3 ) 



for a^a" = 0, 
for CfC 3 > 0. 



(4.24) 



The malonic acid concentration [CH2(COOH)2] is consider "absorbed" into the kinetics of kf B , 

= 2.3 x KT 3 M, 



and the reaction rates are given in Table 4.3 
and the reaction term c 



where [H+] = 0.8 M, [Fe(Phen 



0.43 from (4.9a)— (4.9e). Note that the entropy terms J?^ from §3 present 



a very delicate problem in this nonequilibrium setting. Most explicitly, the equilibrium constant 
itself K eq is not well-defined here, and determining the effective activities fij of the reactions is a 
difficult problem. In our results we merely assume that K eq scales with maxjfc^, which though 
crude, is enough to assure that the entropy conditions from §3 are preserved. Also notice that using 
only a very slightly more complicated approximate form, such as K eq ~ a\j 'a 2_c a 2 a4 immediately 
leads to numerical instabilities due to the concentration scalings in the BZ reaction. 

Let us further comment that it seems that in fact this subtle equilibrium behavior at t n+1 must 
be determined by relying upon a bath assumption at t n that leads to a second order condition on the 
first order assumption employed in the classical Oregonator model. More clearly, in the Oregonator 
model the bath concentrations are assumed large and constant to a first order approximation, 
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Figure 6: On the left is the gradient of a\ (i.e. <Xi) at T = 3 seconds, while the right shows S from 
§3 at T = 3 seconds on Q = [—5, 5] 3 . We emphasize relative magnitudes here, where the gradient 
cr\ goes from highest (red) to lowest (blue), and S is lowest (red) to highest (blue). 

while here, we make no such assumption a priori but must rely upon an equilibrium condition 
that suggests that a\ ^> o^ _c a 2 a4 or at the very least that a\ ^> a^a^, which we consider a 
second order bath assumption because it plays no direct role in the system dynamics, but only on 
the a posteriori entropy consistency and /ip-adaptivity of the solution. 

Here we present some numerical results from the BZ solution, using the initial conditions given 

by: 

C*1,0 = 4 X 10- 5 + 1 X 10-6 e -((-l-6) 2 +( 2 /-1.6) 2 +(^-1.6)^)/2.5 ) 

a 2fi = 1 x 1(T 7 + 1 x i -6 e -((^-i-5) 2 +( 3 /-i.5) 2 +(,-i.5)^)/2. 5) 
a 3 ,o = 2.3 x 1(T 3 + 1 x io-4 e -((-i-4) 2 +(y-i.4) 2 + (,-i.4)^)/2. 5; 
«4,o = 1 x KT 3 , and a 5j0 = 1 x 1(T 3 . 

As we see in Figures [5] and |6j the /^refinement and coarsening is driven by the structure of the 
initial-boundary conditions. We get oscillation behavior of our solutions, and as seen in |4| the 
entropy consistency is preserved and drives the /^refinement regime. 

§4.4 Error behavior at equilibrium 

Finally we consider a simple equilibrium problem comprised of two constituents and constructed 
in such a way as to allow for complete decoupling between the constituents in the mass action, and 
thus obtain an exact analytic solution that may be easily employed for error analysis. We assume 
for this case that 5^ (a) = since the error behavior of the very same LDG method employed here 
has been previously analyzed by the authors in 
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p-convergence in dim=2 



p-convergence in dim=3 
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Figure 7: We plot the p-convergence of the equilibrium solution, where in N = 2 we set h = 1/32 
and for N = 3 we have h = 1/16. 



V 


L 2 -error for iV = 2, h = 32 


L 2 -error for N — 3, h — 16 





0.638527 


2.41962 


1 


0.0230384 


0.283632 


2 


0.000702172 


0.027854 


3 


1.88756 x 10~ 5 


0.00241581 


4 


4.59182 x 10~ 7 


0.000190342 


5 


1.02788 x 10~ 8 


1.38649 x 10~ 5 



Table 4: We give the L 2 -errors shown in Figure ffj 



That is, consider the elementary equilibrium reaction satisfying: 

v{a x v\ot 2 (4.25) 

k b 

such that the coupled system of differential equations is comprised of, 

a[ = v[{k b a U 2 — kfa" 1 ), a' 2 = v\{kja^ — A^c^ 1 ), so that v^a'i = —u(a 2 - (4.26) 

Integrating for any t G [0, T eq ) with T eq the equilibrium time (which exists a priori for min{kb, kf} ^ 
0) and letting the initial concentration a 2 fi = 0, then we further notice that at each t we have 

/ f 

a.\{t) = ol\q \a 2 (t), and at equilibrium that ai(T e „) = a\ o \-a 2 (T eq ). (4.27) 

v 2 ~ v 2 

But then, for the special case of v\ = v{ = 1 assuming ideal behavior where we may take that 
K eq = kf/k b = Oi 2 {T eq ) / ai(T eq ) , where the equations of (4.27) provide that a 2 = cti )0 — act, and 
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/7-convergence in dim=2, p=1 h-convergence in dim=3, p=1 




h level h level 



Figure 8: Here we show the /i-convergence of the equilibrium solution, where p = 1 in both cases. 
For N = 2 the best fit rate of convergence is ~ 2, and in N = 3 is ~ 1.9. 
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p — 1 


p = 2 


L 2 -error 


Convergence Rate 


L 2 -error 


Convergence Rate 


1/4 


1.64322 




0.535138 




1/8 


0.479583 


1.78 


0.0779288 


2.78 


1/16 


0.123667 


1.96 


0.0101416 


2.94 


1/32 


0.0311632 


1.99 


0.00128043 


2.99 


1/64 


0.00780632 


2 


0.000160455 


3 


1/128 


0.00195255 


2 


2.00695 x 10~ 5 


3 


1/256 


0.000488199 


2 


2.50908 x 10~ 6 


3 



Table 5: We give the L 2 -errors and convergence rates shown in Figures |8]-[9] for N = 2. 



h 


p = 1 


p = 2 


L 2 -error 


Convergence Rate 


L 2 -error 


Convergence Rate 


1/4 


3.73802 




0.535138 




1/8 


1.08911 


1.78 


0.0779288 


2.72 


1/16 


0.283632 


1.94 


0.0101416 


2.93 


1/32 


0.0716523 


1.98 


0.00128043 


2.98 


1/64 


0.0179601 


2 


0.000160455 


3 



Table 6: We give the L 2 -errors and convergence rates shown in Figure [8]-[9] for N = 3. 
also yields for the equilibrium constant that K eq = (a 10 — aci(T eq ))/ai(T eq ). Using these relations, 
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Figure 9: The /i-convergence of the equilibrium solution with p = 2. For N 
of convergence is ~ 3, and in N = 3 is ~ 2.9. 



2 the best fit rate 



we 



then rewrite a[ in the first equation of (4.26) as 



ft. 



ki,a2 — hfCti 

hai,o - (h + kf)aci 

(k f + k b )(ai(T eg ) - ax), 



(4.28) 



which has a solution of the form of Q4.7| ) , giving that 

a, = exp" (a 1>0 - ai {T eq )) + a x {T eq ). 



(4.29) 



for any X C [0, T eq ) containing the initial state and any t > T eq , which is just to say the solution 
only depends only upon the initial and equilibrium concentration of a\ — hence fully independent 

of Qf2. 

By contrast we implement our p redic tor multi-corrector in the naive way to recover stf (a) by 



simply solving the discrete form of (4.26) with v\ 
arrive with the solutions 



a 



a 



n+l 



n+1 



exp 



- I At k f dt 



a\ — 



1, such that by the usual procedure we 



K, 



+ 



Oo 



(:(/ 



(4.30) 



exp 



- Iai k b dt 



«2 _ Keq&l) + Keq®l 



where at equilibrium the constant terms balance to unity. 



Then to test our method we compare the error behavior of (4.29) to (4.30) where we take 



an end time for our simulation T which is appropriately set to T > T eq . Here we have a stable 
equilibrium solution (See definition 11.21 in [60]), such that we expect the solution to rapidly 
converge to the equilibrium point in time to machine precision, where the only error remaining 
should be that taken with respect to the standard L 2 -projection. We use the following initial 
conditions: 

a x o = 1 + 4e (*+!) 2 /3.75 ; and a2o = Q ^ 
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All the solutions were run at SSP(5,3) using At = 1 s with T = 200 s. The p-convergence results 
are shown in Figure [7] and Table [4] in both dimension two and three for regular meshes. The h 
convergence results are shown Figures [8] and [9] as well as in Tables [5]-[6] in both dimension two and 
three. Here we use spatially homogeneous refinements of integral value in each direction to obtain 
the expected results. 

§5 Conclusion 

We have developed a predictor multi-corrector time-operator splitting RKLDG SSP scheme that 
utilizes a stability preserving /ip-adaptive entropy consistency scheme for its coarsening and refine- 
ment methodology. The scheme is presented and implemented for arbitrary spatial (i.e. N < 3) 
and component (i.e. n computable) dimension, and includes methods which adopt varying func- 
tional parameters (e.g. *3)(a) and K eq (a)) as well as arbitrary extended Robin boundary data 
as so applied to a generalized subset of react ion- diffusion equations which we denote: quiescent 
reactors. 

The entropy methodology that serves as the fabric of the /ip-adaptive scheme, is extended from 
the regularity analysis of [13], and provides for a sharp stability condition on the computational 
well-posedness of the reaction-diffusion system. 

In addition we have presented solutions to a number of application models, most notably we 
have derived a novel solution to the classical BZ reaction. This is a difficult and complicated 
reaction regime which is often used to underscore the nuance involved in rate-coupled reaction 
mechanisms. 

Our future directions are to extend the quiescent reactors to include fluid reactors where 
density p = p(t, x) and temperature (energy) <£ = <£($) = <£($(£, x)) are fully coupled, in addition 
to adding turbulence models and electromagnetic fields for weakly ionized and plasma reactors. 
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